1 Problema

A motivação original de nosso projeto é o entendimento dos padrões climáticos e suas alterações no clima de uma das cidades mais relevantes num cenário global: Londres. Para isso decidimos fazer um projeto do 1º formato (utilizando um dataset público), onde a base escolhida (london_weather.csv, obtida no website Kaggle através desse link) pode ser interpretada como uma série temporal, a medida que mede as variações do tempo (entre 1979 e 2021) de variáveis fixadas em um único local.

O objetivo do projeto é o uso de um algoritmo de machine learning para que o mesmo consiga prever a categoria de duas das variáveis da database (precipitação e neve) para uma nova informação adicionada ao modelo. As categorias seriam então “chove” ou “não chove” e “neva” ou “não neva”. A ideia será fazer 3 métodos distintos (Regressão Logística, KNN e Árvore de Decisão), compará-los através de suas métricas relacionadas à acurácia e ver qual tem a melhor performance/predição.

2 Dados

A base de dados escolhida foi criada a partir da união de medições oriundas de pedidos de atributos individuais do clima providos pela European Climate Assessment (ECA). As medidas desta base de dados em particular foram gravadas pela estação climática nas redondezas do Aeroporto Heathrow em Londres, Reino Unido. O tamanho original da base de dados escolhida, assim como uma lista dos atributos e suas descrições, está descrito abaixo:

london_weather.csv - 15341 observações x 10 atributos:

  • date - data em que ocorreu a medição - (int)

  • cloud_cover - medição da nebulosidade em oktas - (float)

  • sunshine - medição da luz solar em horas (hrs) - (float)

  • global_radiation - irradiação medida Watt por metro quadrado (W/m2) - (float)

  • max_temp - temperatura máxima registrada em graus Celsius (°C) - (float)

  • mean_temp - temperatura média registrada em graus Celsius (°C) - (float)

  • min_temp - temperatura mínima registrada em graus Celsius (°C) - (float)

  • precipitation - precipitação medida em milímetros (mm) - (float)

  • pressure - pressão medida em Pascals (Pa) - (float)

  • snow_depth - profundidade da neve medida em centímetros (cm) - (float)

2.1 Tratamento dos dados

import pandas as pd
import numpy as np

url = 'https://raw.githubusercontent.com/bmorbin/ML_Project/main/data_raw.csv'
data = pd.read_csv(url)

# Print raw data
print(data.head(5)) 

# Print count of NA for each column
       date  cloud_cover  sunshine  ...  precipitation  pressure  snow_depth
0  19790101          2.0       7.0  ...            0.4  101900.0         9.0
1  19790102          6.0       1.7  ...            0.0  102530.0         8.0
2  19790103          5.0       0.0  ...            0.0  102050.0         4.0
3  19790104          8.0       0.0  ...            0.0  100840.0         2.0
4  19790105          6.0       2.0  ...            0.0  102250.0         1.0

[5 rows x 10 columns]
print(data.isna().sum()) 

# Removing rows with NA, except the column snow_depth
date                   0
cloud_cover           19
sunshine               0
global_radiation      19
max_temp               6
mean_temp             36
min_temp               2
precipitation          6
pressure               4
snow_depth          1441
dtype: int64
data = data.dropna(subset=[col for col in data.columns if col!='snow_depth']) 

# Print count of NA again after delete some NA
print(data.isna().sum()) 

# Applying simple regression model to predict missing values from snow_depth
date                   0
cloud_cover            0
sunshine               0
global_radiation       0
max_temp               0
mean_temp              0
min_temp               0
precipitation          0
pressure               0
snow_depth          1418
dtype: int64
from sklearn.linear_model import LinearRegression 

def reg_simple(data):
    '''
    A function to predict snow depth by mean temperature of the day
    Returns the maximum snow depth value predict
    '''    
    data_reg_train = data.dropna()
    
    X_reg_train = np.array(data_reg_train.mean_temp).reshape(-1,1)
    y_reg_train = np.array(data_reg_train.snow_depth).reshape(-1,1)

    reg = LinearRegression()
    reg.fit(X_reg_train, y_reg_train)
    
    data_NA = data[data['snow_depth'].isna()]
    y_pred = reg.predict(np.array(data_NA.mean_temp).reshape(-1,1))
    return(round(max(y_pred)[0],1))
    
print(f'Minimum snow depth predicted is {reg_simple(data)}')
# So the conclusion is to input 0s in missing values from column snow_depth
Minimum snow depth predicted is 0.2
data.snow_depth[data['snow_depth'].isna()] = 0 

# Print count of NA again after input NA values from snow_depth
print(data.isna().sum()) 

# Convert column date to pandas date type
date                0
cloud_cover         0
sunshine            0
global_radiation    0
max_temp            0
mean_temp           0
min_temp            0
precipitation       0
pressure            0
snow_depth          0
dtype: int64
data['date'] = pd.to_datetime(data['date'].astype(str), format = "%Y%m%d") 

# Creating column month
data['month'] = data['date'].dt.month.astype(str) 

# Creating column rain indicator
data['rain'] = ["1" if p > 0 else "0" for p in data['precipitation']] 

# Creating column snow indicator
data['snow'] = ["1" if p > 0 else "0" for p in data['snow_depth']] 

# Setting the target classifier
data['rain_snow'] = data['rain']+data['snow'] 

# Set column date to index
data = data.set_index('date') 

# Print proportion of time wich rains and snows
print(round(data[['snow','rain']].value_counts(normalize=True),3)) 
snow  rain
0     0       0.515
      1       0.476
1     1       0.005
      0       0.004
dtype: float64
print(round(data.rain_snow.value_counts(normalize=True),3))

# Describe dataframe
00    0.515
10    0.476
11    0.005
01    0.004
Name: rain_snow, dtype: float64
print(data.describe()) 
        cloud_cover      sunshine  ...       pressure    snow_depth
count  15261.000000  15261.000000  ...   15261.000000  15261.000000
mean       5.271018      4.345633  ...  101535.436079      0.034336
std        2.067743      4.019266  ...    1049.536452      0.519855
min        0.000000      0.000000  ...   95960.000000      0.000000
25%        4.000000      0.500000  ...  100920.000000      0.000000
50%        6.000000      3.500000  ...  101620.000000      0.000000
75%        7.000000      7.200000  ...  102240.000000      0.000000
max        9.000000     16.000000  ...  104820.000000     22.000000

[8 rows x 9 columns]

2.2 Dataframe usado

Mostrando as 10000 primeiras linhas do dataframe já tratado.

3 Modelos

Primeiramente devemos encontrar uma solução baseline para o problema proposto. Usando a intuição e o senso comum imagina-se que, em dias muito frios, haverá uma maior probabilidade de nevar, enquanto que dias mais quentes terão maior incidência de chuva e sem neve. Tomemos então isso como nossa solução baseline que será comparada aos modelos categorizados citados acima (Regressão Logística, Árvore de Decisão e KNN).

# 20% for test data
n_test = int(len(data)*0.2)
test_data = data[-n_test:]

# 80% for train and validation
n_train_val = len(data) - int(len(data)*0.2)
df = data[:n_train_val]

target = 'rain_snow'
input_columns = ['cloud_cover','sunshine','global_radiation','max_temp','mean_temp','min_temp','pressure','month']

X, y = df[input_columns], df[target] 

# checking if all classes are present in the training set
print(y.value_counts())
00    6314
10    5764
11      67
01      64
Name: rain_snow, dtype: int64

Nos algoritmos a seguir, serão utilizados 30% do conjunto df (conjunto sem a parte de teste) para validação e 70% para treino.

3.1 Regressão Logística

Para o algoritmo de classificação por Regressão Logística, serão comparados as diferentes combinações de variáveis.

Reference

from sklearn.model_selection import train_test_split
import matplotlib.pyplot as plt
import time
from sklearn.linear_model import LogisticRegression
import itertools

start = time.time()

# creating all combinations of the inputs variables
comb = []
for L in range(1, len(input_columns)+1):
    for subset in itertools.combinations(input_columns, L):
        comb.append(list(subset))
print("Number of combinations: ",len(comb))
Number of combinations:  255
l_total_train = []
l_total_val = []

b=5 #iterations for bootstrap
for c in comb:
    l_a_train = []
    l_a_val = []
    for i in range(1,b):
        X_train, X_val, y_train, y_val = train_test_split(X,y,test_size=0.3, random_state=i)
        X_train_comb, X_val_comb = X_train[c], X_val[c]
                
        clf = LogisticRegression(multi_class='ovr',solver='lbfgs')
        clf = clf.fit(X_train_comb, y_train)
        
        l_a_train += [sum(clf.predict(X_train_comb) == y_train)]
        l_a_val += [sum(clf.predict(X_val_comb) == y_val)]
    
    l_total_train += [np.array(l_a_train)]
    l_total_val += [np.array(l_a_val)]
mean_train = np.transpose(1-np.array(l_total_train)/len(y_train)).mean(0)
sd_train = np.transpose(1-np.array(l_total_train)/len(y_train)).std(0)

mean_val = np.transpose(1-np.array(l_total_val)/len(y_val)).mean(0)
sd_val = np.transpose(1-np.array(l_total_val)/len(y_val)).std(0)

axis_x = list(range(len(comb)))

# reorder the values by the minimum mean error of train to maximum mean error of train
sort_index = [tup[0] for tup in sorted(list(zip(axis_x,mean_train)),key=lambda tup: tup[1], reverse=True)]
mean_train = mean_train[sort_index]
sd_train = sd_train[sort_index]
mean_val = mean_val[sort_index]
sd_val = sd_val[sort_index]

# index of best combination
best_comb_ID = list(mean_val).index(min(mean_val))

# Plotting the combinations with errors
plt.title('Regressão Logística\nErro x Combinação')
plt.xlabel('ID Combinação')
plt.ylabel('Erro') 
plt.plot(axis_x, mean_train, label="Treino")
plt.plot(axis_x, mean_val, label="Validação")
plt.xticks(axis_x, rotation = 90)
plt.fill_between(axis_x, mean_train+sd_train, mean_train-sd_train, alpha=0.5)
plt.fill_between(axis_x, mean_val+sd_val, mean_val-sd_val, alpha=0.5)
plt.axvline(x=best_comb_ID,color='silver', linestyle='dashed')
plt.legend(loc="lower left")
plt.tick_params(axis='x', which='major', labelsize=1)
plt.show()

end = time.time()

print("Spending time to Logistic Regression analysis of combinations:",(end-start)/60, "minutes")
Spending time to Logistic Regression analysis of combinations: 3.7618120551109313 minutes

3.2 Árvore de Decisão

Para o algoritmo de classificação por Árvore de Decisão, comparou-se os valores de profundidade máxima da árvore que apresentaram menor erro de validação.

Reference 

from sklearn import tree

l_total_train = []
l_total_val = []

start = time.time()

b=20 # iterator bootstrap (resampling the set)
c=20 # number of complexity based on max_width
for j in range(0,b):
    X_train, X_val, y_train, y_val = train_test_split(X,y,test_size=0.3, random_state = j)
    
    l_a_train = []
    l_a_val = []
    
    for i in range(1,c):
        clf = tree.DecisionTreeClassifier(max_depth = i, criterion = "gini")
        clf = clf.fit(X_train, y_train)
        l_a_train += [sum(clf.predict(X_train) == y_train)]
        l_a_val += [sum(clf.predict(X_val) == y_val)]
    
    l_total_train += [np.array(l_a_train)]
    l_total_val += [np.array(l_a_val)]

mean_train = np.transpose(1-np.array(l_total_train)/len(y_train)).mean(1)
sd_train = np.transpose(1-np.array(l_total_train)/len(y_train)).std(1)

mean_val = np.transpose(1-np.array(l_total_val)/len(y_val)).mean(1)
sd_val = np.transpose(1-np.array(l_total_val)/len(y_val)).std(1)

# best max_depth argument to don't overfit
best_depth = list(mean_val).index(min(mean_val))+1

axis_x = list(range(1,c))

# Plotting 
plt.title('Árvore de Decisão\nErro x Complexidade')
plt.xlabel('Complexidade: Profundidade Máxima')
plt.ylabel('Erro') 
plt.plot(axis_x, mean_train, label="Treino")
plt.plot(axis_x, mean_val, label="Validação")
plt.xticks(axis_x)
plt.fill_between(axis_x, mean_train+sd_train, mean_train-sd_train, alpha=0.5)
plt.fill_between(axis_x, mean_val+sd_val, mean_val-sd_val, alpha=0.5)
plt.axvline(x=best_depth,color='gray', linestyle='dashed')
plt.legend(loc="lower left")
plt.show()

end = time.time()
# time spending to load all algorithm analysis
print("Spending time for Decision Tree analysis of complexity:",(end-start)/60, "minutes")
Spending time for Decision Tree analysis of complexity: 0.43555124203364054 minutes

3.3 KNN ( K-Nearest Neighbors )

Notebook Reference

from sklearn.neighbors import KNeighborsClassifier
b=20 # iterator bootstrap (resampling the set)
c=20 # number of complexity based on number of neighbors
for j in range(0,b):
    X_train, X_val, y_train, y_val = train_test_split(X,y,test_size=0.3, random_state = j)
    
    l_a_train = []
    l_a_val = []
    
    for i in range(1,c):
        clf = KNeighborsClassifier(n_neighbors=c)        
        clf = clf.fit(X_train, y_train)
        l_a_train += [sum(clf.predict(X_train) == y_train)]
        l_a_val += [sum(clf.predict(X_val) == y_val)]
    
    l_total_train += [np.array(l_a_train)]
    l_total_val += [np.array(l_a_val)]

mean_train = np.transpose(1-np.array(l_total_train)/len(y_train)).mean(1)
sd_train = np.transpose(1-np.array(l_total_train)/len(y_train)).std(1)

mean_val = np.transpose(1-np.array(l_total_val)/len(y_val)).mean(1)
sd_val = np.transpose(1-np.array(l_total_val)/len(y_val)).std(1)

best_nn = list(mean_val).index(min(mean_val))+1

axis_x = list(range(1,c))

# Plotting for view the best number of neighbors
plt.title('KNN\nErro x Complexidade')
plt.xlabel('Complexidade: número de vizinhos')
plt.ylabel('Erro') 
plt.plot(axis_x, mean_train, label="Treino")
plt.plot(axis_x, mean_val, label="Validação")
plt.xticks(axis_x)
plt.fill_between(axis_x, mean_train+sd_train, mean_train-sd_train, alpha=0.5)
plt.fill_between(axis_x, mean_val+sd_val, mean_val-sd_val, alpha=0.5)
plt.axvline(x=best_nn,color='silver', linestyle='dashed')
plt.legend(loc="lower left")
plt.show()

end = time.time()

print("Spending time for KNN algorithm analysis:",(end-start)/60, "minutes")
Spending time for KNN algorithm analysis: 3.720773720741272 minutes

4 Comparando modelos

A partir da seleção de parâmetros de cada algoritmo realizada nas seções anteriores, compara-se então os modelos entre si para analisar qual possui melhor acurácia.

from sklearn.metrics import plot_confusion_matrix

X_train, X_val, y_train, y_val = train_test_split(X,y,test_size=0.3, random_state = 0)
X_train_comb, X_val_comb = X_train[comb[best_comb_ID]], X_val[comb[best_comb_ID]]

print("Accuracy of Regression Model: ",sum(LogisticRegression(multi_class='ovr',solver='lbfgs').fit(X_train_comb, y_train).predict(X_val_comb) == y_val)/len(y_val))
Accuracy of Regression Model:  0.6033306033306033
plot_confusion_matrix(LogisticRegression(multi_class='ovr',solver='lbfgs').fit(X_train_comb, y_train), X_val_comb, y_val)  
<sklearn.metrics._plot.confusion_matrix.ConfusionMatrixDisplay object at 0x0000000049EB7EB0>
plt.show()

print("Accuracy of Decision Tree: ",sum(tree.DecisionTreeClassifier(max_depth=best_depth, criterion = "gini").fit(X_train, y_train).predict(X_val) == y_val)/len(y_val))
Accuracy of Decision Tree:  0.7447447447447447
plot_confusion_matrix(tree.DecisionTreeClassifier(max_depth=best_depth, criterion = "gini").fit(X_train, y_train), X_val, y_val)  
<sklearn.metrics._plot.confusion_matrix.ConfusionMatrixDisplay object at 0x000000004AA3CFD0>
plt.show()

print("Accuracy of KNN: ",sum(KNeighborsClassifier(n_neighbors=best_nn).fit(X_train, y_train).predict(X_val) == y_val)/len(y_val))
Accuracy of KNN:  0.6743106743106743
plot_confusion_matrix(KNeighborsClassifier(n_neighbors=best_nn).fit(X_train, y_train), X_val, y_val)  
<sklearn.metrics._plot.confusion_matrix.ConfusionMatrixDisplay object at 0x0000000005CDB280>
plt.show()

5 Conclusão

Conclui-se então que o algoritmo árvore de decisão apresenta maior acurácia, além de ser o algoritmo mais rápido para análise. Portanto, roda-se a seguir o modelo com o conjunto de dados de teste:

# Fit the model selected
model = tree.DecisionTreeClassifier(max_depth=best_depth, criterion = "gini").fit(X_train, y_train)

X_test = test_data[input_columns]
y_test = test_data[target]

print("Accuracy of Decision Tree to test data: ",sum(model.predict(X_test) == y_test)/len(y_test))
Accuracy of Decision Tree to test data:  0.7260812581913499
plot_confusion_matrix(model, X_test, y_test)  
<sklearn.metrics._plot.confusion_matrix.ConfusionMatrixDisplay object at 0x000000004C7362B0>
plt.show()

Todos códigos referentes ao projeto podem ser encontrados nesse repositório.